Field theoretic calculation of the surface tension for a model electrolyte system 
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' We carry out the calculation of the surface tension for a model electrolyte to first order in a 

' cumulant expansion about a free field theory equivalent to the Debye-Hiickel approximation. In 

CNj , contrast with previous calculations, the surface tension is calculated directly without recourse to 

integrating thermodynamic relations. The system considered is a monovalent electrolyte with a 
region at the interface, of width h, from which the ionic species are excluded. In the case where 
, the external dielectric constant eo is smaller than the electrolyte solution's dielectric constant e we 

■ show that the calculation at this order can be fully regularized. In the case where h is taken to be 
' zero the Onsager-Samaras limiting law for the excess surface tension of dilute electrolyte solutions 

is recovered, with corrections coming from a non-zero value of eo/e. 

o 

; I. INTRODUCTION 

, . The first experiments to measure the surface tension of electrolyte solutions show that the excess surface tension, 

■ denoted in this paper by (Je, due to the presence of the electrolyte is positive 0. This result has been confirmed 
Q , by more recent experiments 0- This effect was explained by Wagner Q who pointed out that when the dielectric 
Q ' constant of the bulk solvent (here water) e is greater than that of the exterior (here air) eo then the image charges, 

due to the dielectric variation across the surface, repel the solute ions from the surface and thus lead to a reduction 
of the density of ions near the surface with respect to the bulk. Applying the Gibbs adsorption isotherm we then find 
that de must be positive. In addition experimental results on systems at weak dilution for solutes of the same valency 
I are very similar, suggesting a universal limiting law at weak dilution. Such a universal limiting law was subsequently 
\^ ' obtained by Onsager and Samaras 

A series of experiments carried out in the 1930s caused a certain controversy as at very small electrolyte 
concentrations a negative excess surface tension was reported. It seems that these experiments have not been revisited 
using modern techniques, or at least have not been reproduced since. If a negative Ue is found, then appealing to the 
^0 Gibbs adsorption isotherm, there must be some mechanism causing ions to be positively adsorbed near the interface. 
_ Various authors have discussed ion-specific effects which could explain such a phenomenon |E 0, S S ^| i and also 
. lead to the ion dependent variations seen in the measurements of cTg at higher concentrations. 
' The calculation of the surface tension of electrolytes was recently revisited in a series of papers by Levin 
^ , and Levin and Flores-Mena Because of the thermodynamic equivalence of ensembles, an exact calculation of the 
'■^ ' surface tension should give the same result independent of the ensemble chosen since the thermodynamic identities from 
^ ' which the surface tension is calculated are exact. However, Levin points out that calculations of the surface tension 
O ' invariably rely on approximation schemes, notably the Debye-Hiickel approximation, and that a given approximation 
[ scheme will generally yield different results for different choices of thermodynamic ensemble. For example. Levin 
^ . applies a canonical approach whereas the original Onsager-Samaras result was obtained using the grand canonical 
ensemble. In the approach of Levin de is given by the excess Helmholtz free energy due to the presence of an interface. 
This free energy excess is obtained by calculating the internal energy due to the presence of the interface and then 
integrating it via the Giintelberg charging process to obtain the free energy. In the limit of weak electrolytes the 
Onsager-Samaras limiting law is recovered thus, as Levin remarks, suggesting that the Onsager-Samaras limiting law 
is indeed exact. 

In this paper we calculate de in the grand-canonical ensemble by directly calculating the excess grand potential 
due to the presence of an interface. In this way we avoid the integration of differential thermodynamic identities 
such as the Gibbs adsorption isotherm or the Giintelberg charging process, and so provide another route for doing 
the calculation. In addition, we develop a controlled per turbation theory based on a cumulant expansion, similar to 
that used for bulk electrolytes by Netz and Orland [1^; this is a perturbation expansion in the coupling constant 
g = Ib/Id, where Id is the Debye length and Ib the Bjerrum length. We show that the Onsager-Samaras limiting law 
is the first term in this cumulant expansion, showing that it is indeed exact to this order. The limiting laws obtained 
in the literature are given in the limit where eo/e ^ 0, which is clearly a good approximation for aqueous solutions 
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in air where eo/e ~ 1/80. In this paper we generaUze the Onsager-Samaras result and give the corresponding hmiting 
law in the case where eg/e > 0. 

Our approach is also applied to a modified model of the interface where there is surface-exclusion layer for the 
ions of thickness h: a region at the surface from which the hydrated ions are forbidden 

II in El- Highly accurate 

numerical integration is used to investigate the importance of the effect of this exclusion layer on the value of Ce- 

The techniques used in this paper are based on the field theoretic Sine-Gordon representation of the grand partition 
function first introduced in this context in [Tsj l . The perturbation theory about the free field or Debye-Hiickel theory 
is carried out using a functional path integral technique introduced recently by the authors 16, 17, 18], which lends 
itself to the geometry of planar systems and gives a powerful alternative method for the calculation of the functional 
determinants involved. 

We conclude with a discussion of our results and the possible advantages of our approach for calculating in more 
complex models where, for example, a surface charge exists due to a thermodynamic adsorption process for one of 
the ionic species at the surface or due to a difference in the hydrated radii between the cations and anions. 



II. THE MODEL 



We consider a model consisting of a semi-infinite electrolyte bulk with monovalent salt in contact with a semi-infinite 
exterior, see Fig. The bulk solvent's dielectric constant is denoted by e and the exterior dielectric constant is 

denoted by sq. There is a region of width h between the exterior and bulk which is filled with the bulk solvent but 
from where the salt ions are excluded, this is a standard surface-exclusion layer and was first introduced by Randies 
|l4| in the context of electrolyte surface tensions. The width h of the surface-exclusion layer is the order of a hydrated 
ion radius. For simplicity both hydrated anions and cations are taken to be of the same size and are hence both 
excluded from this region and so there is no surface charging process. However, the model and approach can be 
generalized to ions of different radii which will lead to different ion-specific surface layer widths and so allowing a 
charging mechanism |17| . In the bulk solution the fugacity of the cations and anions is equal and denoted by fi. 
The system up can be summarized in terms of a spatially dependent dielectric constant e(z) and spatially dependent 
fugacity fi{z) which are defined as follows 

e(z) — So z < —h 

e{z) = e z > -h (1) 

and 

fi{z) ^ z <0 

H{z) ^ n z>0 (2) 



In the grand canonical ensemble the grand partition function for the system is given by the functional integral over 
the Wick rotated electrostatic potential (p 



S = y d[0]exp(5[</>]) (3) 

with the action S given by 

S[(t)] = y dx e(x)(V0)2 + 2jdx ^(x)cos (e^0) (4) 

where e is the electron charge and /3 is the inverse temperature. We take the area of the system across the interface 

to be A and the length (in the z direction) of the exterior to be L and of the bulk to be L' . If one considers just the 
exterior system without any interface, its grand partition function is given by 

^E = J dMexp(5£[0]) (5) 

with Se given by 

Se[cP] = ^^ I d^eo{Vq>f (6) 
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FIG. 1: Schematic image of the exterior bulk interface for the model considered here. The values of the local dielectric constants 
and fugacities as a function of the distance from the dividing surface are shown. The charges of the ions are taken to be at 
their centers which are excluded from the surface-exclusion layer of width h. 

the integration being over the region L' x A. For a pure bulk system with no interface the grand partition function 
is given by 



5s - J d[cl)]exp{SBm 



where Sb is given by 



dx eiycjjf + 2 dx ^cos (e/3(/)) 



(7) 



(8) 



the integration being over the region L x A. The surface tension is then given by the difference in the grand potential 
of the system with the interface and that of the sum of the two individual (exterior and bulk) systems divided by the 
total area, i.e. 



a=^{J{L',L)^j(^\L)-j(^\L')) 



(9) 



where J^^^ = — ln(5£;)//3 and J^^^ = — ln(SB)//3 denote the grand potentials for a bulk system of electrolyte and 
exterior system of the same volumes [L x A] and [L' x A] respectively, but with no interfaces. The definition of Eq. 
for the surface tension is of course also in agreement with various other methods for calculation. For example it 
is the same as that obtained from the Gibbs adsorption equation as originally used by Onsager and Samaras. The 
expression Eq. Q for the surface tension can also be obtained from the formula 



2(7 = - 



Pd{L)dL 



(10) 



where Pd{L) is the disjoining pressure for a film of external medium of thickness L surrounded by bulk electrolyte |l9j . 
This system consists of two bulk surfaces a distance L apart, and so twice the surface tension is given by the work 
needed to create an infinitely thick film: L — > cx). As mentioned in the introduction, the approach here is different 
from previous techniques since here the grand potential difference corresponding to the surface tension is calculated 
directly. 

The excess surface tension ae for a system with bulk electrolyte concentration p is defined by 



<7e{p) = a{p) - a(0) 



(11) 
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where (t(0) is the surface tension of the system with no added electrolyte. This definition means that (Jg is free of the 
ultraviolet or short-distance divergences found in calculations of the surface tension between two media of differing 
dielectric constants |3l23- 

In electrostatic problems where the chemical potential and dielectric constants depend only on the coordinate z, the 
field theory can be formulated as a functional path integral for a dynamical field 4>{r, z) which evolves in a temporal 
coordinate z |18| |. The functional Hamiltonians are denoted by He in the exterior region, Hb in the bulk and Hs in 
the surface-exclusion layer. In three dimensions this functional problem cannot be solved exactly but in one dimension 
it can be and leads to an explicit solution for the one-dimensional Coulomb gas ^1^. The free Debye-Hiickel theory 
can be also solved in this formulation and one can develop a perturbation theory about it as we shall show here. 
For the moment we will use the Hamiltonian formulation explicitly in order to find a formal expression for the excess 
surface tension. For a globally electro-neutral system with no interfaces and Hamiltonian H and length L in the z 
direction, one may write the grand partition function as 

S = Trexp(-LiJ) (12) 

that is, we take the system to be periodic in the z direction. Hence for the pure bulk of electrolyte density p one has 
that for large L 

S(^) = {^i''\p)\exp{-LHB{p))\^lf\p)) (13) 

and for the exterior region 

S(^) = (vI/(^Vxp(-i'i?i^)|*D (14) 

where |*[,^^^''^) and l^^f^) are the normalized ground-state wave functionals for the bulk and exterior functional 
Hamiltonians Hb (p) and He respectively. Note that the wave functionals must be normalized so that the correspond- 
ing grand potential is zero for a system of zero volume, that is, zero length in the z-direction. If the corresponding 

( B) ( E) 

ground-state energies are ' {p) and , then we have 

/3J(^)(L) ^ l4^'(p) (15) 
PJ^^\L') = L'Ef^ (16) 



and the corresponding bulk pressures are given by 



/3P(^) = (18) 

For the system with interface we find 

S = (vI/(^) I exp {-L'He) exp {~hHs) exp (-(L - h)HB)) \^f\p)) (19) 
this is easily seen by joining two such systems together with periodic boundary conditions. We thus obtain 

/3J(L,L') ^ + {L- h)E^^\p) In ((vj/f ^ | exp (-/li/^) \^[''\p)))] (20) 

Using this the excess surface tension is given by 



/^(i?^(p)-4''^(0))+ln 



^[''\e^v{-hHs)\^f\p)) 



(vI.(^)|exp(-/.iJs)|*r(0)), 



(21) 



Using the relations Eq. (|17|l and Eq. H18|l we thus obtain 
where 

AP(p)=Pb(p)-Pb(0) (23) 

is the bulk pressure due to the presence of the electrolyte. The expression Eq. H22I) is difficult to evaluate, although 
an approach using standard quantum-mechanical perturbation theory might be investigated. However, if the original 
field theory is free or Gaussian, Eq. (|22|l is relatively straightforward to compute. We shall use Eq. (|22|) to evaluate 
the contribution to the surface tension coming from the cumulant expansion about the free Debye-Hiickel theory. 
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III. CUMULANT EXPANSION OF THE EXCESS SURFACE TENSION 



Perturbation theory about the Debye-Hiickel theory is carried out by decomposing the action S in the following 
manner 



S = So + AS 



(24) 



The first term 5*0 is a Gaussian or free term given by 



So = - — 



f3e 



l-L'-h,-h]xA ^ 
dx ((V<?i)2 +771^02) 



dx (V0)2 



-h,()]xA 

2fi{L-h)A 



l[0,L-h]xA 

where m is the Debye mass given by = 2pe'^(3/e. The correction to the Gaussian action AS* is given by 



AS-: 



[0,L-h]xA 



2fi (cos (e/3(/)) - 1) 



(25) 



(26) 



The term AS* is of order of the dimensionless coupling constant g = Ib/Id where Id ~ \/ ^/m is the Debye length and 
Ib = e'^ 13/ Aire is the Bjerrum length. A cumulant expansion in AS" generates a resummed expansion in g in the sense 
that the term of order n in the cumulant expansion has the form C„ = g"fn{g)- In the bulk the function fn{g) then 
has the form fn{g) = X]m=i O'n.mg"^- However in the presence of the interface we will see that fnig) has an extra 



term containing logarithmic terms in g of type ^ 



m=l "'n,m 



(7™ln(g). This can be shown by considering the form of 



the bulk action Sb written in times of the dimensionless field 0' — eP/cft^/g and by measuring length in units of the 
Debye length (y — mx). In the new field and length variables one has the bulk action 



l\2 



dy ^(Vc^') 



zjg) 

Ang 



dy cos (V50') 



(27) 



where Z{g) is given by Z{g) = 1/ (cos {^(f)')). Here we have used the fact that at a point x in the system, the average 
density of cations / anions is given by 



/0±(x) = A*(exp (±ie/?(/)(x))) 

and for x in the bulk jO±(x) = p. It is easy to check that Z{g) = 1 + zig + Z2g^ 
the bulk as above we obtain 



Sb — So 



AS 



(28) 

Using the same decomposition in 



(29) 



where 



is the Gaussian or free action and 



AS^ 



Airg 



dy 



Z{g)cos {^gct^') + U'^ - Z{g) 



(30) 



(31) 



Using the series form for Z{g) we see that AS" can be expressed as a power series in g with first term 0{g). It can 
also be shown that (AS') = at 0{g) for the homogeneous bulk system; the corollary is that (AS") ^ at 0{g) only 
for systems which are not translationally invariant such as the system with an interface under discussion here. The 
outcome is that when calculating to leading order in g we just need to keep the first term in the cumulant expansion 
of the free field theory with AS treated as a perturbation and the 0{g) contributions to J'^^\L) and j'^-'(L') in Eq. 
are zero. We write 



S = y d[(f\ exp (5*0 + AS) « exp {{AS)o) J d[^] exp {So) 



(32) 
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with 



(A5)o 



/ d[(j)]AS expjSo) 
Jd[cl,]exp{So) 



(33) 



This first term in the cumulant expansion can also be shown to begin with the two-loop term of the standard loop 
expansion, and hence we shall also refer the calculation that follows as the two-loop, or more correctly, the resummed 
two-loop calculation. 

To this order of approximation the grand potential is given by 



J = Jo + A J 



with 



-13 Jo = ln(^y exp (5o) 
-/3AJ = {AS)o 



(34) 

(35) 
(36) 



The action Sq is Gaussian and we define the correlation function of the field at the same point and a distance z 
from the surface-exclusion layer, by 



(0(r,z)0(r,z)>o =G(0,z) 



(37) 



where we have used the fact that the system is isotropic in the plane A {r Cz A) but is not isotropic in the direction 
z. As the action 5*0 is purely quadratic we also have that 



{cl>{r, z))o = 



(38) 



For the bulk system (i.e., without an interface) we note that the same-point field correlator at this level of approxi- 
mation is given by 



(0(r,z)0(r,z))o = Gb{0) = G(0,oo) 



(39) 



since the physics as z — s- oo for the system with an interface at z = is the same as that of the bulk system. Using 
this result, we find that for the system with interface 



/•oo 

(AS')o = A / dz 
Jo 



^ , . e2/32G(0,z), \ ^e^i'^.n ^' 
2m exp( ' ) - 1 + ^^G(0, z) 



Using Eq. (|28|l to relate p and /i, we find that to 0{g) the fugacity is determined by 



p ~ H exp 



-G(0,oo; 



Using the results above, we find 
{AS)o = Ajdz 
where we have defined 



„ / , e2^2^(0,z), , e^(3^G{0,oo) 

2p exp( — 



-) - exp(- 



Gr{0,z) = G(0,z) - G(0,c») 



G(0,z) 



(40) 



(41) 



(42) 



(43) 



Since we seek a result accurate to 0{g) we may expand the second exponential in the integral in Eq. H42II to first 
order and neglect higher 0{g^) terms. Using the definition of the Debye mass m this yields 



{AS), 
A 



J dz 2p ^cxp(— 



e2/32Gfl(0,z) 



Gfl(0,z) 



The first term in Eq. (|44ll is finite even in the limit h ^ {) whereas the second term 



r = 



dzGR{0,z) 



(44) 



(45) 
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is ultra-violet divergent as h 0. This divergence is due to the integral over the potential due to the image charge. 
We might naively resolve this potential difficulty by observing that if we also expand the first exponential in Eq. H44|l 
this term is exactly cancelled (AS')o = to 0{g), so resolving the difhculty. However, this is expansion is incorrect 
since this divergence is, in fact, cancelled by another arising in Jq. The expansion of the first exponential gives rise 
to an erroneous divergence which then survives wrongly in the final result; there is no such divergence. The form of 
Eq. (|44(l is familiar since the first exponential is the Boltzmann factor for the repulsive image charge potential that 
we should expect to appear and is reminiscent of terms in the the Mayer expansion. 

To calculate G(0,z) it is convenient to use the path integral representation of the problem. Using 

<P{i^,z) ^ -^^4>{P,z)exp{ip-r) (46) 
^ p 

we find that the Gaussian action 5*0 simply becomes sum of independent Harmonic oscillators 

So = 2fiL + Y^ Sp (47) 



where 

(48) 



Sp — —— I dz 



M(.)^^ + M(«V'(P,«W(P)«-P) 



where M{z) = f3e{z), uj{p,z) = \p\ — p for z £ [—L',h] and uj{p,z) = -y/p^ + m^ for z e [h,L]. By expanding in 
terms of the Fourier modes we find that 

^(0, z) - jY.(Hp, z)4>hp, z))o (49) 
p 

The Euclidean Fe ynm an propagator for a simple Harmonic oscillator, with Hamiltonian denoted by Ho{uj, M), over 
a time t given by |22 | 

{X\exp(-tHo(uj,M))\Y) = f ^{^ , x 1 exp {--Mlo coiUuot) [X^ + _ 2Xrsech(wi)] ) (50) 

yzTTsmh [ujt) J \ 2 J 

and the ground-state wave function is given by 

{M^,M)\X) = (^^^ ' exp (^-h'luX^^ (51) 

with energy Eo{oj, M) = lli/2. In the free field theory we thus find 

,7. x^ _ iM^Ejp), Me)\ exp i-hHojLJsip), Ms)) exp {-zHo{ujb{p), Mb)) X^\M^b{p), Mb)) , . 

mP,zm P.z))o {^,(,,^(p),ME)\exp{-hH,iLUsip),Ms))expi-zH,iLOBip),MB))\M^B{p),MB)) ^ ' 

where the subscripts B, E, and S refer to the bulk exterior and surface-exclusion layer values of the various simple 
harmonic oscillator Hamiltonians Ho and the corresponding masses M and frequencies w in these regions. 
Carrying out the Gaussian integrations we thus obtain that 

(^(p,z)0(-p,2))o = i?33' (53) 

where D is the matrix 



a 


-b 





-b 


c 


-d 





-d 


e 



D ^ \ -b c -d\ (54) 



The elements of D are given by 

a — [3eop + (Sep coi\i{ph) 

b — fiep cosech{ph) 

c = (3ep coth(pft,) -I- [3t\J p^ + coih{\J p^ + m? z) 



d = lit\/p'^ + m^cosech(\/p^ -I- vn? z) 

e = I3e\/ p'^ + m?{l + coth(\/p2 -|- rn? z)) (55) 
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A long but straightforward calculation now gives the result 



G(0,z) 



dk k 



Kcoth{Kmz) + kB 
K{kB + K){1 + coth{Kmz)) 



(56) 



where the integral over k is between and A/m where A is an ultra-violet cutoff in the Fourier modes of the field 
in the plane A. In the present calculation we will see there are no ultra-violet divergences and we may take the limit 
A ^ oo. In Eq. (|^ and through out the rest of this paper we use the following definitions 



and 



where 



B = 



K = ^JW+l 

A cxp(— 2fcm/i) 



1 



1 + Acxp(-2fcm/i) 



A = 



Using Eq. we find that 

and using Eq. and Eq. we obtain 



G(0,cx)) 



m 
47r/3e 



dk 



K 



„ k{K~kB) , , 

dk—, V exp(~2Kmz) 

K{kB + K) ' 



A{zm) 



(57) 



(58) 



(59) 



(60) 



(61) 



In the case A = 1, Levin and Flores-Mena in Eq. (8) of reference ^3 quote a similar formula for W(z) in their 
notation. Comparing our result at A = 1 with theirs, we note a misprint where the exponential exp(— 2fc(z — d)) in 
the integrand of their equation should read exp(— 2pz). With this correction we identify 



W{z) = |A(mz) 



with d = h. 



A=i 



(62) 



Our result, however, applies for all A, < A < 1 and all h > 0. 
Using Eq. and Eq. (gSJ we find 



P9_ 
2m 



dk k 



K -kB 
K^{kB + K) 



(63) 



Repeating the above calculation for a pure bulk system, we see that in the absence of an interface that Gii{z, 0) = 
and consequently that the corresponding term (AS')o is zero, and so for the pure bulk without interface we have to 
one loop that J^^^ — Jq^^ ■ For a pure exterior system the action is purely Gaussian and AS**-^^ — identically, and 
so to one-loop Eq. Q becomes 



The excess surface tension is thus given by 



-(AJ + Jo- J(f-^(f) 



ae(p)=a:(p)+a(")(p) 



(64) 



(65) 



where 



(66) 
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and (7e°^ (p) is the excess surface tension for a system with just the action 5*0 which can be calculated exactly in the 
quantum mechanical formulation as all the simple harmonic oscillators are decoupled. We have from Eq. (|22|) 



af\p) = 2^^/!--^^ [/i(i?o(c^i3(p,p),MB)-i?o(coB(p,0),MB)) 
p 

/ (V^o(^i;(p), Me)\ exp {-hHojusiv), Ms)) \M^b{v, P),Mb)) \ 
. {M^e{v),Me)\ exp {~hHo{uos{v), Ms)) \M^b{v, 0), Mb)) ) 



(67) 



where we have made explicit the dependence of the bulk frequencies lob on p, wb(p;p) = ^Jp^ + m?{p). Note that 
the first term in the right-hand side of the above comes from the constant term in the action 5*0. 
Using Eqs. H50|l and (|51|l we obtain 



af\p)^PDebyeh+-^ / kdk 



K — k \ K 

21n ( 1 + ^, (1 + Aexp(-2fcmft.)) - ln(— ) 
2k Ik 



(68) 



with Poebye the Dcbyc pressure, that is to say the bulk pressure to 0(.g), given by 



2p 



47r 



kdk{K - k) 
^2p(l-|) 



2p 



247r 



where the rightmost expression in Eq. H69|l is obtained after calculating /i in terms of p |l6j| . 
Collecting all these contributions we arrive at our final result for the excess surface tension 



2ph [1 



2p 



d{mz) 



1 — exp 



-A(mz) 



^2m ' '^^ ^ 



, , K-k, , , , ,,A , (kB-K) 
41n (1 + ^(1 + AcM-2kmh)) j - 21n(-) + ^ j,^ 



(69) 



(70) 



where the function A{mz) as defined by Eq. (j61|l . and we have arranged the terms to explicitly show the dependence 
on the dimensionless coupling g. we denote the first term to be the exclusion term, the second to be the depletion 
term and the third to be the Casimir term. 

To evaluate this expression it is convenient to decompose A{mz) into a component which is singular as z — ^ 0, 
which gives the direct interaction with the image charge, and a component finite in this limit: 



A{mz) 



Aexp(-2m(z + h)) 
2m(z + h) 

de sinh e 



/ exp(-2mzcosh 9 - 29) (1 - exp(-4m/i sinh 9)) 
V 1 + A exp{-2mh sinh 9 - 29) 



A exp(— 2to2 cosh 9) {exp{—2mh sinh 6) — exp{—2mh sinh 9)) 



(71) 



where the change of variable k — cosh 9 has been used. 

The result for Ce is correct in perturbation theory to 0{g) and holds for < A < 1 and h > 0. In the depletion 
term, the function A{mz) is the potential due to the interaction of a charge with its image and, as is seen above, not 
only includes the screened Coulomb (Yukawa) potential, which is singular as z — > when /i = 0, but also contains 
non-singular correction terms which, in particular, are important when h > 0. We have derived (|70|) directly from the 
perturbation expansion for the free energy but the same result would be obtained from the Gibbs adsorption isotherm 
or the Giintelberg charging process; in both cases a perturbation expansion can be obtained for the appropriate 
quantity which is then appropriately integrated. Levin and collaborators have derived a similar result to Eq. 
(|7n)l at 0{g) for the case A = l,h = but they assume the phenomenological form for A{niz) given by the screened 
Coulomb potential at ft. = 0: the first term in Eq. 1)7111 . As we shall see in the next section, the result by Levin 
for cTg is numerically similar to ours when evaluated at A = l,/i = but for general values of A, ft the full result 
for A{mz) in Eq. (|71fl is needed for an accurate calculation of the depletion term. The Casimir term is generated 
automatically in the Giintelberg charging process used by Levin but again to obtain the general result correct to 0{g) 
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presented here, the process must be derived from the perturbation expansion for the energy density considered as a 
function of the electric charge e. In addition, in our approach, whatever the method for deriving Ue, the perturbation 
series for ae can systematically calculated to higher orders in g by including terms of higher order using the cumulant 
expansion in AS", Eq. (|^ . 

In the next two subsections we discuss the consequences of this result. 



A. The Onsager-Samaras limiting law 



In this section we shall consider the case where /i = and the case e > eg. We show how the Onsager-Samaras 
limiting law 4] for (Tg at A = 1 follows from our result and we derive the generalization to cases where A < 1. 
When /i = Eq. ifTHjl becomes 



2p 



d{u) 



1 — exp 



-A{u) 



— I 

4m 



(72) 



We thus see that the calculation of cre{p) to the first order in the cumulant expansion about the Debye-Hiickel ap- 
proximation is divergence free. We now discuss the physical origins of these terms. The first term gives a contribution 
to the surface tension ai^'^ which can be interpreted as being proportional to the depletion of solvent with respect to 
the bulk at the interface within the Debye-Hiickel approximation. This term appears in the original Onsager-Samaras 
calculation where it is then integrated with respect to the fugacity via the Gibbs adsorption equation to obtain the 
excess surface tension. In the Debye-Hiickel approximation it is easy to see that 



(3a, 



{D) 



dz (p+(z) - 2p) , 



(73) 



where p±{z) indicates the average cation/anion density at a distance z from the surface. 
When ft, = we have 



^(y) ^ ^^^P( + (1 _ a2) / de sinh e exp {-2u cosh i 
2w Jo 

We find the asymptotic expansion of in the limit of small g to be 



exp(-26') 
1 -f Aexp(-26 



(74) 



P9^ 
2m 



1 



2A: 



■(1 + A) (2Aln(2)- (l-fA)ln(l-f A)) 



0(9^ Hg)) 



(75) 



When A = 1 Eq. (|75|l is in agreement with the result of Onsager and Samaras 3 , thus showing that the limiting 
law is exact up to the order of the correction indicated in Eq. (|75|l . We note that from our earlier discussion higher 
order corrections coming from the cumulant expansion will also be 0{g^). 



B. The general case 

Our results for o-g and A(mz) in Eqs. H70|) and (|71|l apply generally for all < A < 1, /i > 0. When h is non-zero 
the addition of another length scale in the problem renders the derivation of analytical results considerably more 
complicated. The first term of Eq. H7()(l has a simple physical interpretation, it gives a contribution Poebyeh to ae 
which can be interpreted as the work done to expel the ions from the surface-exclusion layer into the bulk. In the 
limit where hm <C 1 «. e. h <^Id the second two terms of Eq. (|70|l we can set ft « and recover the Eq. (|75|l for 
these two terms. We now present some numerical results based on the highly accurate Vegas |23| integration package. 

To carry out the integration over z in Eq. H70|l we need to accurately determine the function A(rnz) given in Eq. 
(|61|l and this itself requires an integration over k. The decomposition for A{mz) given in Eq. (|71|l is vital for good 
convergence of this latter integral since it is converted to a well behaved integral over 6 with the singular nature 
of the potential A{mz) expressed explicitly. To attempt the integration over k in Eq. (|61|l numerically would not 
accurately produce this singular behaviour especially in the region where it is most important, namely as z — > 0. We 
use Vegas to carry out the dO integration in Eq. H71|) and so accurately determine A{mz) on a discrete set of closely 
spaced points for z in the range < mz < 4 and use interpolation to evaluate this function at intermediate points. 
To calculate ae we carry out the separate integrations in Eq. (|70|l again using Vegas. Accurate convergence of the 
numerical integration in all cases is rapid and errors are negligible. In Fig. |(2J) we show the separate contributions 
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to (Te of the exclusion term, the depletion term and the Casimir term (respectively the first, second and third terms 
in Eq. JTU))) and the total value, as a function of solute molarity < a; < 1.0 for h = 0.0, 0.1, 0.2, 0.3 (nm) and a 
temperature of 20°C and A = 0.975, appropriate for water where e/eo ~ 80. 

From Fig. ||2J) we note that for < h < 0.3(nm) the dependence of cTg on h is rather mild especially at low solute 
density, for example 0.2(mol); ae initially decreases with increasing h but then increases as the exclusion contribution 
begins to dominate and the effect of the depletion and Casimir terms is reduced. Obviously, for larger h the domination 
of the exclusion term is complete and ae with rise linearly with h, Eq. (|70|) . for fixed solute density. However, the 
range of h considered here is typical of physical films and should be compared with the Debye length at solute density 
of l(mol) and T = 20°C of Id = 0.305(nm). 

It is interesting to compare the result for the depletion term calculated from the full expression for A(mz) given 
in Eq. I|71|l with that for A{mz) approximated by the first term: the Yukawa potential. This is relevant because 
the Yukawa contribution has an obvious physical significance as the potential for the image charge repulsion and 
is the extension to non-zero h of the potential used by Levin ('|ll|'). We can then examine the importance of the 
non-singular correction term in Eq. (|71|l . whose origin is no so phenomenologically obvious. For solute density of 
l(mol) we show in Table HJ the respective contributions of these two calculations as h increases for two values of 
A = 0.975, 0.6 at T = 20°C. The first value corresponds to the water-air interface with e/Zao/eo = 80, f-exterior/ (-o = 1, 
and the second value is for an interface between water and an exterior medium with ^exteriorl^o — 20. In each case, 
the first column gives the contribution when A{mz) is is approximated by the Yukawa term and the second column 
tabulates the contribution when the full expression is used for A{mz). For the water/air interface there is negligible 
difference for h = 0(nm) but whilst both contributions decrease with h the full result is over five times larger than the 
phenomenological Yukawa approximation suggests. The difference is much more marked in the case with A = 0.6 even 
at h = 0(nm) with the full result an order of magnitude larger than the Yukawa approximation when h ~ 0.3(nm). 
Note that the Debye length is Id = 0.305 (nm), comparable with the largest value of h here. These results show that 
in a realistic film, which will generally have a surface layer of thickness in the range discussed here, the corrections to 
the Yukawa approximation to the image-charge interaction Eq. I|7HI are overwhelmingly important, especially when 
the thickness h > Id- For higher solute densities this inequality likely to be easily satisfied. A similar effect occurs 
for the Casimir term, and it is the slower decrease of both these terms with h compared with the phenomenological 
prediction that almost exactly balances the increase of the exclusion term with h so that the h dependence of ae is 
relatively weak in the range shown in Fig. (for A = 0.975 here). An outcome is that the result of Levin for 
(Te, which applies only to the case h = 0(nm) A = 1, is numerically similar to ours for A = 0.975 but here we have 
extended the results accurately to general A and h to 0{g). 

It should be noted that in the exclusion term in Eq. 170|l the Debye-Hiickel formula for the pressure has been 
used and for solute density of l(mol) and T — 20°C the dimensionless coupling constant is g — 2.33913, and so 
the Debye-Hiickel correction to the free gas law pressure is nearly 40%. This indicates that corrections to the bulk 
pressure at 0{g^) and higher will make a significant contribution at this and higher solute densities and, by inference, 
the higher order corrections to both the depletion and Casimir terms in Eq. (|70|l should be calculated. This is the 
aim of work in hand. 

In Fig. (PJ we show the temperature dependence 10°C < T < 30°C for different values of h and solute density of 
l(mol). It is clear from the h — 0(nm) results that the temperature dependence of the depletion and Casimir terms 
is very weak, and that the dominant contribution for h > 0(nm) is from the exclusion term and simply comes from 
the T-dependence of the free gas pressure, (x T, plus the dependence of the Debye-Hiickel correction cx T~^/^. 

IV. CONCLUSIONS 

We have shown that the calculation of ae for a simple model of an electrolyte can be formulated in terms of a 
perturbation expansion in the dimensionless coupling constant g = Ib/Id, where Ib and Id are the Bjerrum and 
Debye lengths, respectively. For the first time we derive the full general expressions for ae to 0{g) for general values 
of A = (e — eo)/(e + ^o), < A < 1 and h > 0(nm). The calculational method is based on a direct calculation of 
the grand-potential difference between a system with bulk/exterior interface and a bulk and exterior system with no 
interface (both with periodic boundary conditions). In this simple model an exclusion layer for the hydrated ions 
at the surface was included, both cations and anions we implicitly taken to be of the same size and thus had the 
same range of exclusion h. Due to the symmetry between cations and anions in the model here, no mean-field or 
average electrostatic potential or effective surface was generated. The Onsager-Samaras limiting law is shown to be 
the limiting form of the first term in this expansion for small g and for the first time we have derived its generalization 
in Eq. H75I) to the case where A < 1, h = 0(nm). 

The method is equivalent to other approaches to calculating CTe such as the Gibbs adsorption isotherm and the 
Giintelberg charging process, both of which can be formulated using our techniques as perturbation series in g. The 
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FIG. 2: The component contributions to the surface tension in mN/m versus the molar density for surface exclusion layer 
thickness h = 0.0, 0.1, 0.2, 0.3(nm), A = 0.975 and T = 20°C. 



strength of our method is that it gives a systematic expansion which can be extended to higher orders in g by including 
the higher-order terms in the cumulant expansion in AS in Eq. H31|l . The expUcit terms given here in Eq. (|70|l to 
0{g) are respectively the exclusion term due to the exclusion of the gas of solute ions from the surface layer, the 
depletion term which gives the contribution from the image-charge repulsion for ions approaching the surface, and 
the Casimir term arising from the change in energy of electric field modes due to the presence of the surface. Terms 
higher order in g will correct the first two of these contributions. Our method also gives the exact form for these 
terms, and especially gives the full expression for the image charge potential G(z, 0) = q A(mz)/e'^0 defined in Eqs. 
(|61|l and (|71|l . The expected screened Coulomb (Yukawa) potential used by Levin '11] can be identified from Eq. 
(|71|l but the non-singular correction term is not so easily argued phenomenologically, and from Table (P) is seen to be 
important for /i > 0. By inference a similar effect occurs for the Casimir term. For the exclusion term in Eq. H70|l . 
dominant for large h, higher order corrections to the Debye-Hiickel approximation are necessary for solute densities 
greater tha n l( mol), and by inference the other terms will also need correction at higher densities. It should be noted 
that Levin [ll| uses the free gas law to compute the exclusion term omitting the Debye-Hiickel correction which is 
equivalent to setting = in this term. 

The formalism used here can be extended to deal with cases where an effective surface charge is present due to 
either a difference in hydrated ionic sizes; the behavior of the surface charge in such a model was recently analyzed in 
a weak charging approximation by the authors |17|. In addition one may also apply this formalism to other systems 
with different energetic or thermodynamic mechanisms leading to surface charging; these are the so-called charge 
regulated models 24]. In all cases, one can go beyond the first order expansion used here, although this will require 
a more sophisticated theory with a short distance cut-off to regularize the divergences arising at higher orders in 
perturbation theory; for example, the use of a regularization scheme based on an additional repulsive short range 
Yukawa interaction will permit the use of the path integral techniques employed in this work. 
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A = .975 


A= .6 


h{nm) 


Yukawa (mN/m) 


Full (mN/m) 


Yukawa (mN/m) 


Full (mN/m) 


0.0 


0.710 


0.712 


0.561 


0.598 


0.1 


0.286 


0.414 


0.188 


0.331 


0.2 


0.106 


0.287 


0.067 


0.241 


0.3 


0.042 


0.231 


0.026 


0.205 



TABLE I: The contribution to ae of the second (depletion) term in Eq. I70II for solute density of 1 (mol) as a function of 
the exclusion layer thickness h in (nm). The second column gives the contribution when A{mz) is defined in Eq. 17H is 
approximated by the first (Yukawa) term and the third column tabulates the contribution when the full expression is used for 
A{mz). The results shown are for two values of A = 0.975, 0.6 and T — 20°C. The first value corresponds to the water-air 
interface with eH2o/eo = 80, texterior/^o ~ 1, and the second value is for an interface between water and an exterior medium 
with eexterior/^o ~ 20. For the water/air interface there is negligible difference for h — 0(nm) but whilst both contributions 
decrease with h the full result is over five times larger than the phenomenological Yukawa approximation suggests. The 
difference is much more marked in the case with A = 0.6 with the full result being an order of magnitude larger than the 
Yukawa approximation when h — 0.3(nm). Note that the Debye length is Id ~ 0.305(nm), comparable with the largest value 
of h here. 
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FIG. 3: The temperature dependence of a for solution of density one mol for T in range 10 — 30(°)C and for surface exclusion 
layer thicknesses h = 0.0, 0.1, 0.2, 0.3(nm). The dependence on T becomes more marked as h increases, being about 10% over 
this range for h = 0.3(nm). This effect is almost entirely due to the exclusion contribution whose T-dependence comes from 
the formula for the free gas pressure oc T plus the dependence of the Debye-Hiickel correction oc T~^^^ . Here A = 0.975. 
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